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ABSTRACT 

Many questions in physical cosmology regarding the thermal and ionization history of the intergalactic 
medium are now successfully studied with the help of cosmological hydrodynamical simulations. Here 
we present a numerical method that solves the radiative transfer around point sources within a three 
dimensional cartesian grid. The method is energy conserving independently of resolution: this ensures 
the correct propagation speeds of ionization fronts. We describe the details of the algorithm, and compute 
as first numerical application the ionized region surrounding a mini-quasar in a cosmological density field 
at z = 7. 

Subject headings: radiative transfer — galaxies: formation — cosmology: theory — galaxies: 
intergalactic medium — galaxies: ISM 



1. INTRODUCTION 

The incorporation of radiation processes and radiative 
transfer into cosmological hydrodynamical simulations is 
essential for modeling the structure and evolution of the Ly 
a forest (Cen et al. 1994; Zhang, Anninos & Norman 1995; 
Hcrnquist et al. 1996; Hachnclt & Stcinmetz 1997), inter- 
preting the observations of helium and metal line absorbers 
at high redshift (Rcimers et al. 1997; Hellsten et al. 1997; 
Songaila 1998), and for simulating the early reheating and 
reionization of the intergalactic medium (IGM) (Ostriker 
& Gnedin 1996; Gnedin & Ostriker 1997). In the last 
case the large computational effort is justified by the hope 
that the study of the transition from a neutral IGM to one 
that is almost fully ionized could provide some hints on the 
first generation of stars and quasars in the universe. Un- 
til now, radiative transfer effects have either been ignored 
in such simulations, or treated as a self-shielding correc- 
tion to an optically thin approximation (Katz et al. 1996; 
Gnedin & Ostriker 1997). Typically, radiation fields are 
treated as isotropic backgrounds J v {z) which are either 
specified as an external function computed by other means 
(e.g., Haardt & Madau 1996), or are computed by aver- 
aging over all sources within the computational volume 
(Cen 1994; Gnedin & Ostriker 1997). Only approximate 
treatments exist of the reprocessing of radiation via ab- 
sorptions internal to the sources as well as by the IGM, an 
effect that can significantly influence the spectrum of the 
metagalactic flux (Haardt & Madau 1996). Moreover, at 
z > 3 (5) the IGM itself is opaque in the helium (hydrogen) 
Lyman continuum, and radiation backgrounds become 
increasingly inhomogeneous and anisotropic (Reimers et 
al. 1997; Madau, Haardt & Rees 1998). Recent ana- 
lytic work has attempted to forge a closer connection be- 
tween sources, transport, and sinks of cosmic radiation, 
albeit in a spatial- and angle-averaged way (Haiman & 
Loeb, 1998a, b; Madau, Haardt & Rees 1998). 



In this Letter we describe an approach to cosmological 
radiative transfer which relaxes these assumptions, and is 
appropriate on scales comparable to the separation of indi- 
vidual sources of radiation. Our method has the property 
that energy is conserved independently of numerical reso- 
lution, ensuring, for example, that ionization fronts prop- 
agate at the correct speeds. The radiation driven front 
surrounding a mini-quasar in a cosmological density field 
at z = 7 is computed as a first application of the technique. 

2. COSMOLOGICAL RADIATIVE TRANSFER 

The equation of cosmological radiative transfer in co- 
moving coordinates (cosmological, not fluid) is (Norman, 
Paschos & Abel 1998): 



c dt 



n ■ VI„ 



H{t) dl v 
c dv 



%I v) =t]v- Xvlv (1) 



where I v = I(t, x, h, v) is the monochromatic specific in- 
tensity of the radiation field, h is a unit vector along the 
direction of propagation of the ray; H(t) = a/ a is the 
(time-dependent) Hubble constant, and a = is the 

ratio of cosmic scale factors between photon emission at 
frequency v and the present time t. Here r\ v , \vi and c 
denote the emission coefficient, the absorption coefficient, 
and the speed of light, respectively. Equation 1 will be 
recognized as the standard equation of radiative transfer 
with two modifications: the denominator a in the second 
term, which accounts for the changes in path length along 
the ray due to cosmic expansion, and the third term, which 
accounts for cosmological redshift and dilution. In princi- 
ple, one could solve equation (1) directly for the intensity 
at every point in (t, x, h, v) space given r\ and \- However 
the high dimensionality of the problem not to mention the 
high spatial and angular resolution needed in cosmological 
simulations make this approach impractical. Therefore we 
proceed through a sequence of well-motivated approxima- 
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tions which reduce the complexity to a tractable level. 

Comparing the third with the second term in equation 
(1), the classical transfer equation (Kirchhoff 1860): 



1 dl v 

— + n ■ VI„ =T) V - Xulu 
c at 



(2) 



is recovered if the scale of interest, L, is much smaller 
than the horizon, c/H(t). (Note that this is only true if 
I 1 ' 7^7 < Iu\ for continuum radiation. However, we 
can nonetheless still use eq. (2) for line transfer provided 
we Doppler shift the absorption cross section in \v due 
to Hubble expansion.) In the case of constant absorption 
and emission coefficients the time derivative can also be 
dropped and one is left with the static transfer equation 



h ■ VI V = r\ v - Xulv 



(3) 



This equation is also adequate for problems in which the 
absorption and emission coefficients change on timescales 
much longer than the light crossing time, L/c. This will al- 
ways be the case in the volumes we will be able to simulate 
in the near future, and thus we adopt this approximation. 
At small distances from the source, however, this simplifi- 
cation will always break down allowing I-fronts to expand 
faster than the speed of light. This can be easily dealt 
with as we will show in the subsequent section. 

An extensive literature exists on methods of solution 
for equation (3) (e.g., Mihalas & Mihalas 1984) depending 
on the symmetry of the problem and the properties of r\ v 
and Xv- A direct solution of equation (3) is impractical 
due to the high dimensionality of the problem. Our ap- 
proach is approximate, and is based on decomposing the 
radiation field into point source and diffuse components: 
I v = If, ts + . The motivation for doing this is that 
whereas If, ts is highly anisotropic, we expect to be 

nearly isotropic, and accordingly we can employ the best 
methods taylored for each component. 

Utilizing the linearity of the radiation field, we can 
write: 



(4) 
(5) 



where diffuse emission (e.g., due to recombination radia- 
tion) appears only in the transport equation for the dif- 
fuse radiation field via the source function S v = f] v /xv 
This has the advantage that equation (5) can now be 
solved ray by ray in a completely decoupled fashion. For 
Xv{x) = const, equation (5) has the simple analytic solu- 
tion 

I v {x\) = I„{x )e~ T = I„(x ) exp(-Xv(xi ~ x )). (6) 

In the next section we describe an efficient implementa- 
tion for solving equation (5) on rays emanating from point 
sources within a uniform Cartesian grid. Several meth- 
ods of solution for equation (5) are discussed in Norman, 
Paschos & Abel (1998). 

3. IMPLEMENTATION 

Figure 1 gives the flowchart of the simple algorithm; de- 
tails of the individual computational steps arc given in the 
following sections. For each point source in our volume, a 



set of radial rays quasi-uniformly distributed in solid angle 
are constructed. Enough rays are used such that every cell 
at large distances from the source is crossed by at least one 
ray on average. We are able to use fewer angles than one 
would naively assume by Monte Carlo sampling in angle 
within a radiation-matter coupling timescale. Each ray is 
discretized "cast" into ray segments according to how it 
intersects the cell boundaries. On each ray segment the 
intensity is attenuated according to equation (5) using an 
absorption coefficient appropriate to that cell. The num- 
ber of absorptions in the frequency interval v-\-dv in each 
ray segment inside a cell is simply related to the decrease 
in I'P ts along that segment. The total number of photoion- 
izations in the cell is the sum over all ray segments crossing 
the cell. 



Fig. 1. — The flowchart for the radiation transfer algorithm. 

3.1. Choosing Angles 

To properly describe the I-front one needs at least 
N a = f a x 2irr/ Ax rays so that at least one ray is cast to 
each cell of sidelength Ax at the equator of a sphere with 
radius r. We introduce the factor f a to allow us to control 
the number of rays. To adopt the minimum number of 
rays needed we store the maximum radius, r max , needed 
to capture the furthest point of the I-front. 

We use spherical coordinates (r, <j>, 6) , 



x = r cos 4> cos 9, < (j> < 2ir 
y = r sin <j> cos 9 
z = r sin 9, 



7T 7T 

- <9< - 
2 ~ ~ 2 



(7) 



and divide the sphere in segments of similar size such that 
all rays will have similar fluxes. We can now find the an- 
gles that define the rays which pierce through segments of 
roughly equal area. The discrete values for 9 are given by 



1 x 2tt tt N a 

9j = (j--) , 1 < 7 < — 

3 u 2' N a 2 ~ J ~ 2 



(8) 



Closer to the poles fewer azimuthal angles are required, 
Nl = max(N a cos 0,1). (9) 



The azimuthal angles are now chosen to be, 



(10) 



One also needs to know the area of the sphere segment, 
A(i,j) described by each ray. In units of the surface of the 
sphere (4wr 2 ), this is 



Mh 3) 



sin 92 — sin 9\ | 
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(11) 
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Fig. 2. — Sphere from top with the division into areas of similar 
size. The azimuthal angle < <j> < 2-k and 59 are also indicated. 
The variying N$(6) is evident. 

A random rotation of the coordinate system is then in- 
troduced for each snapshot solution of I v to avoid any 
possible pole artifacts and reduce the number of angles 
required. 

3.2. Casting Rays 

Ray casting is the next step. This means to compute, 
for a given ray (0, <fi), the indices k(l) of the grid 

cells it traverses, as well as the lengths of the ray seg- 
ments dS(l) within each cell. Here k are the indices 
of grid cells in a Cartesian lattice x(i),y(j), and z(k), and 
/ is the index of the ray segment. For simplicity we re- 
strict our discussion to uniform, isotropic grids such that 
x(i) = iA, y(j) = jA, z{k) = kA where A is the cell size 
and < i < Ni, and similarly for j and k. Ray casting is a 
problem in computational geometry with efficient methods 
discussed in books on computer graphics. For a uniform 
Cartesian grid the problem is straightforward. Here we 
write down a direct algorithm (i.e., requiring no costly 
sorts). Consider a point source at coordinate x s ,y s ,z s in 
cell i s ,j s ,k s . We have simply = i s ,j(l) = j s ,k(l) = 
k s . We cast a ray in direction 9,(f> and ask which cell 
boundary it intersects first. If it intersects an x boundary 
first, i(2) = i(l)+sign(cos</>), j(2) = j(l),fc(2) = fc(l), and 
analogously if the y or z boundaries are crossed first. To 
determine which cell face is crossed first, we compute the 
distances (radii) from the source to the next x, y or z cross- 
ings. Call these distances r x ,r y and r z . The minimum of 
these three distances determines which cell boundary is 
crossed first. 

This can be expressed as follows: Define 

fi = cos 4>, s M = sign(^i) 
7 = sin0, s 7 = sign( 7 ) 
C=sin0, Sf = sign(C) (12) 

then, for all 2 > 2 we have 



r B (I) = |A[i(0 + -(s„ + 1)] - x s |/max(e, |/icos0|) 
r v (l) = \A[j{l) + -(s 7 + 1)] - t/ s |/max(e, | 7 cos0|) 
r z {l) = \A[k{l) + i(a c + 1)] - z s |/max(e, | CI) (13) 
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S(l) = r y (l) 
j(l)=j(l-l) + s 7 
if MO < min(r x Q),r y (l))], then 
S(l) = r z (l) 

k{l) = k(l - 1) + s c . (14) 
Here e is a small number to avoid dividing by zero. 

3.3. Computing Rates 

Although the above technique can be applied to any ab- 
sorption process, in the following we will discuss the spe- 
cific case of hydrogen photoionization, and compute the 
photoionization rate coefficient and the associated heating 
term. Additional species can be treated in an analogous 
way. Given the path length dS(l) = S(l) — S(l — 1), and 
the absorption coefficient in cell (i(l),j(l),k(l)), the op- 
tical depth for photoionization along this path is given 
by T i — Xu{i(l),j{l),k{l))dS(l). Further, from the so- 
lution of the static transfer equation (6), one sees that 

St x N„(l — 1)(1 — e~ r ") photons are absorbed in the time 
interval, St. Hence, the rate of change of the neutral hy- 
drogen density, nn, due to photoionization is simply given 

by 



n H = N l ~ 



1 - e 



/Vc 



cell, 



(15) 



where N u is the monochromatic number of photons emit- 
ted by the source per unit time, and Vceii denotes the vol- 
ume of the cell. The heating rate is then, 



e ={hv-l ryd) N y 



l-e 



/V c 



cell- 



(16) 



Both of the above rates are summed over the contributions 
of all rays. Employing the analytical solution for each ray 
segment insures that I-fronts move at the correct speed 
independent of spatial resolution. 

3.4. Faster than light I-fronts 

In deriving equation (3) we have neglected the time de- 
pendent term of the transfer equation (1). The assump- 
tion that the light crossing time is shorter than the ion- 
ization timescale breaks down close to the source. As a 
consequence the I-front expands at a speed exceeding the 
speed of light. This can also be seen from the simple jump 
condition in a static universe, 



2 drj 

AnrrUii = N — Airas 

dt 



n e n p r 2 dr, (17) 



where as is the recombination coefficient to the excited 
states of hydrogen. In this expression, dri /dt exceeds the 
speed of light for radii less than r c — (N /Aircn-ft) - 5 « 

5.3 kpc(A56/nH) ' 5 - To avoid this unphysical effect we 
simply do not compute rates at distances further than 
ri = c(t — ton). The radius n is also used in the com- 
putation of the optical depths in equations (16) and (15). 
Additionally, to speed up the calculation for radii r < r c 
we evolve the equations on a tenth of the light travel time 
across one cell. As a consequence the time evolution of 
the ionized fration close to the source will not be com- 
puted correctly. However, we do not consider this a severe 
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limitation since the neutral fraction becomes practically 
zero in regions close to the source within a light travel 
time. 

4. A TEST AND AN APPLICATION 

In this section we first provide a test case that demon- 
strates the accuracy of the method, and then, as first appli- 
cation, compute the ionized region surrounding a "mini- 
quasar" in a cosmological density field at z = 7. Mini- 
quasars could be common at early epochs if black holes 
can form in the first 10 s M Q cold dark matter (CDM) 
condensations (Haiman & Loeb 1998b). 

4.1. Spherical I-front Test 

To test our algorithm we set up a point source in a uni- 
form medium, and neglect the effect of radiative recombi- 
nations. The radius of the spherically expanding I-front 
can be derived analytically by balancing the number of 
emitted photons with the hydrogen atoms that are present 
in the ionized volume, 




timels'l 



Fig. 3. — Comparison of the analytical radius of an I— front [eq. 
(18)] with the numerical solution on a 64 3 3D cartesian grid. The 
error bars indicate the maximum deviation found on the spherical 
ionization front. The deviations are always less than one grid cell. 
The straight line indicates were part of the I-front leaves the simu- 
lation volume. The cell size is 1 kpc in this test case. 

Figure 3 compares this analytical solution with the results 
of our algorithm on a 64 3 grid. The radius and the spher- 
ical geometry are perfectly recovered within the accuracy 
of the spatial resolution. The specific parameters used 
in this test were, N = 10 51 s _1 , rin = 10 -2 cm -3 , and 
Ax = 1 kpc. The algorithm has also been tested to repro- 
duce accurately the size of classical Stromgren sphere in 
calculations where radiative recombinations are included. 

4.2. Cosmological Density Field 

Our first interesting application is the evolution of the 
ionization zone surrounding a mini-quasar which turns 
on at z = 7 in the inhomogcneous density field derived 
from cosmological simulations. For simplicity, we have ne- 
glected any thermal or dynamical evolution of gas, and 
included hydrogen radiative recombinations by assuming 
a constant Case B recombination rate of as = 3.4 x 1CP 13 
cm 3 s _1 everywhere on the grid. As a background medium 
we use the H I distribution computed in a SCDM cosmol- 
ogy (with n b = 0.06 and h = 0.5) at z = 6.941 from 
Bryan et al. 1998. The 2.4 Mpc box length corresponds 
to 300 proper kiloparsecs. At the densest cell, which is 
found in a virialized halo of total mass w 1.3 x 10 11 M Q , 



we introduce a quasar-type source with an ionizing pho- 
ton emission rate of N — 5 x 10 53 s _1 . The gas clumping 
factor is («-h)/(^h) 2 = 59.2, and is larger than the ionized 
hydrogen clumping factor, C = («Hii)/( n Hii) 2 , as clumps 
that are dense and thick enough to be self-shielded from 
UV radiation remain neutral and do not contribute to the 
recombination rate. 




Fig. 4. — Propagation of an R-type I— front in a 128 3 cosmological 
density field produced by a mini-quasar with N = 5 X 10 53 s — 1 . The 
solid contours give the position of the I-front at 0.15, 0.25, 0.38, and 
0.57 Myr after the quasar has switched on. The underlying greyscale 
image indicates the initial H I density field. 

Figure 4 illustrates the evolution of the I-front during the 
first 0.6 Myr. The initial spherical expansion (at the speed 
of light since the number of available photons exceeds by 
far the number of neutral hydrogen atoms and recombina- 
tions) is quickly broken as the front quickly expands first 
into the voids (increasing their thermal pressure by many 
orders of magnitude by photo-heating), then more slowly 
into the denser filaments. 

In a highly inhomogcneous universe, the volume- 
averaged gas recombination timescale, 

tree = (UeaBCy 1 = 0.1 GjX C 20 (19) 

(for Ofe/ilo = 0.06) is much shorter than the then Hubble 
time. At later epochs, when the size of the H n region is 
large compared to the scale of the clumping, the front will 
fill its time-varying Stromgren volume in a few recombi- 
nation times just like in the static case (Madau, Haardt, 
& Rees 1998), 

Vs = ^ = 063Mpc3 ^ 53 Ji + A 6 c -i (20) 

n H V ° / 

This is 20 times larger than our simulation box. 

5. SUMMARY 

Not only for theoretical reasons, but also in light of 
upcoming space missions such as MAP, PLANCK, and 
NGST, a detailed understanding of the thermal history 
of the IGM and reionization is highly desirable. In this 
paper we have described a photon conserving algorithm 
which allows to simulate inhomogenous reionization self- 
consistently within in a cosmological hydrodynamical sim- 
ulation. The method employs an on the spot approxima- 
tion to treat the effects of diffuse emission of ionizing pho- 
tons. For scenarios in which stellar sources are responsible 
for hydrogen reionization the on the spot approach is ex- 
pected to give reliable results. However, for calculations 
in which diffuse radiation due to recombination of helium 
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is also expected to be important a separate solver can 
be used. In applications where radiative recombinations 
can be neglected the number of rays used can be reduced 
greately due to the variable choice of coordinate system 
for the selection of rays. Recently, a differnt algorithm for 
3D radiative transfer in cosmological situation has been 
presented by Razoumov and Scott (1998). Their explicit 
advection scheme (at the speed of light) seems well suited 
for situations in which the I-fronts move faster than any 



hydrodynamical flow. Only temporary 2D arrays are used 
in the method presented here requiring negligible small ad- 
ditional memory for 3D cosmological hydrodynamics sim- 
ulations. 
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